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Abstract 


An  approximate  solution  is  presented  to  the  seismic  inverse 
problem  for  two-dimensional  velocity  variations.  The  solution  is 
given  as  a multiple  integral  over  the  data  observed  at  the  upper 
surface.  An  acoustic  model  is  used  and  the  reflections  are  assumed 
to  be  sufficiently  weak  to  allow  a "linearization"  procedure  in  the 
otherwise  non-linear  inverse  problem.  Synthetic  examples  are 
presented  demonstrating  accuracy  of  the  method  with  dipping  planes 
at  angles  up  to  45°  and  with  velocity  variations  up  to  20%.  The 
method  was  also  tested  under  automatic  gain  control,  in  which  case 
velocity  estimates  were  lost  but  the  method  nonetheless  sucessfully 
migrated  the  data. 
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Introduction 


In  recent  papers  (Bleistein  and  Cohen,  1977;  Cohen  and  Bleistein, 
1977),  the  authors  have  shown  that  closed  form  approximate  solutions 
for  the  velocity  (or  a velocity  related  quantity)  profile  can  be 
obtained  for  a wide  variety  of  wave  propagation  equations.  "Closed 
form,"  as  used  here,  means  that  the  velocity  profile  (depth  section) 
is  obtained  by  direct  processing  on  the  observed  data  (time  section) 
itself;  i.e.,  by  performing  weighted  quadratures  on  the  data. 

In  this  paper  we  present  an  approximate  solution  to  an  inverse 

problem  often  used  in  seismic  modeling.  It  is  assumed  that  the 

velocity  in  the  subsurface  varies  in  only  two  dimensions,  vertically 

and  laterally.  A line  of  point  sources  is  set  off  and  the  backscattered 
+ 

signal  is  observed.  This  is  the  mathematical  idealization  of 
summation  of  common  depth-point  (CDP)  traces.  An  approximate  integral 
equation  is  derived  for  the  variation  in  velocity  from  some  known 
reference  value.  When  that  reference  value  is  a constant,  the 
integral  equation  has  a closed  form  solution  for  the  two-dimensional 
velocity  variation. 

The  simplifying  assumptions  which  lead  to  the  closed  form  solu- 
tion presented  here  are  as  follows.  Firstly,  it  is  assumed  that  the 
acoustic  wave  equation  adequately  describes  the  subsurface  wave  pro- 
pagation (although  in  the  first  paper  cited  above  inversion  for 

t 

With  a planar  array  of  backscatter  observations,  the  analysis  can 
be  extended  to  yield  an  approximate  solution  to  the  three-dimensional 
inverse  problem;  i.e.,  to  find  a velocity  that  varies  in  all  three 
dimensions. 

1 


•t  • 


•♦.If-*  '*■  f 


elastic  wave  propagation  was  discussed).  Secondly,  the  subsurface 
velocity  variations  must  be  small.  This  justifies  linearization  of 
an  otherwise  non-linear  problem.  Even  this  "smallness"  limitation  is 
not  overly  restrictive  as  can  be  seen  from  one  of  the  examples  below 
in  which  a 20%  velocity  variation  was  successfully  migrated  and 
estimated.  In  fact,  the  real  world  data  restrictions--noise,  attenu- 
ation, discretization  and  finiteness  of  observations,  etc., —are 
usually  of  greater  concern  than  this  theoretical  linearization 
assumption. 

This  second  assumption  has  signigicant  theoretical  consequences, 
namely,  that  reflecting  interfaces  cause  "weak"  reflections  and  the 
(downward)  transmitted  wave  may  be  taken  to  be  the  response  to  the 
source  in  a homogeneous  (constant  velocity)  medium. 

A computer  program  was  developed  to  implement  our  result  on 
synthetically  generated  data.  The  synthesis  embodies  several  real 
world  restrictions,  namely: 

(i)  the  observations  are  made  only  at  discrete  points 
on  the  line; 

(ii)  the  observations  are  made  only  over  a finite  length 
of  the  line; 

(iii)  the  observations  are  band  limited  in  frequency. 

For  the  examples  here,  in  (i)  the  spacing  between  shot-receiver 
points  was  taken  to  be 

AC  * 100  ft.;  (1) 
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in  (ii),  the  line  was  taken  to  extend  from  -L  to  L,  with 

L=  4,000  ft.;  (2) 

in  (iii)  the  bandwidth  f_  f was  used  with 

= 4f_  = 24  hz.  (3) 

An  important  feature  of  the  direct  inverse  procedure  used  here 
is  that  it  does  not  break  down  as  the  dip  angle  of  reflecting  sur- 
faces is  increased.  This  can  be  seen  in  the  examples  below  of 
dipping  planes  with  dip  angles  up  to  45°.  This  is  so  because  this 
direct  inversion  procedure  is  formulated  from  propagation  governed 
by  the  wave  equation  itself,  rather  than  from  a parabolic  approxi- 
mation to  the  wave  equation.  Parabolic  approximations  are  known  to 
"flatten"  the  dip  of  wave  fronts  (Claerbout,  1976)  and  consequently 
flatten  the  dip  of  profiles,  as  well.  This  is  demonstrated  with 
three-dimensional  examples  in  the  above-cited  paper  by  Bleistein  and 
Cohen  (1977)  in  which  a three-dimensional  inversion  procedure,  starting 
from  the  parabolic  approximation,  is  presented. 

This  feature  of  accuracy  at  all  dip  angles  is  also  present  for 
migration  schemes  based  on  the  wave  equation  itself--see,  for  example, 
French  (1975),  Lamer  and  Hatton  (1976),  Schneider  (1970).  However, 
as  migration  techniques,  these  methods  produce  only  a subsurface 
profile— i .e. , layer  mapping  but  no  velocity  estimates.  The  velocity 


profile  produced  by  the  direct  inversion  method  described  here,  pro- 
duces velocity  estimates  as  well  as  a layer  mapping.  These  velocity 
estimates  require  true  amplitude  information.  When  such  information 
is  not  avai lable--e.g. , if  automatic  gain  control  has  been  applied-- 
the  direct  inversion  technique  simply  migrates  the  data  to  produce  a 
depth  section  without  a velocity  estimate. 


The  linearization  procedure  used  in  the  derivation  of  our 
inversion  procedure  is  often  referred  to  as  the  Bom  approximation 
(Morse  and  Feshbach,  1953).  It  has  been  successfully  exploited  in 
the  past  to  yield  approximate  solutions  of  other  inverse  scattering 
problems.  For  example,  Prosser  (1969,  1976)  has  shown  that  in 
certain  cases,  this  approximation  leads  to  a first  iterate  in  an 
iteration  scheme  to  solve  certain  "refraction"  or  "potential" 
inverse  problems.  Three  essential  requirements  of  his  proof  are 
that  (i)  the  scattering  potential  or  reflectivity  be  "sufficiently 
weak;"  (ii)  the  scattering  region  be  finite  in  extent  and  (iii)  the 
scattering  potential  or  tractive  index  be  "smooth"  in  a manner 
prescribed  by  certain  decay  estimates  on  its  spatial  Fourier  transform 
at  infinity. 


-• 


The  condition  (i)  is  assumed  by  us,  as  well.  Conditions  (ii) 
and  (iii)  are  simply  not  true  in  seismic  problems.  In  particular, 
condition  (ii)  leads  to  a particularly  simple  first  approximation  of 
the  index  of  refraction,  namely,  that  its  three-dimensional  spatial 
Fourier  transform  is  proportional  to  the  backscattered  data.  The 
lack  of  transverse  confinement  heads  to  the  much  more  difficult 
inversion  formula  given  by  equation  (9),  below. 
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Analysis 


It  is  assumed  that  the  medium  to  be  probed  supports  acoustic 
waves  with  wave  speed  v(x,  z).  Here  x is  the  transverse  variable 
and  z is  the  depth  variable  (measured  positive  downward  from  the 
upper  surface).  Thus,  the  medium  has  velocity  variation  in  one 
transverse  direction  only. 

The  governing  equation  for  the  wave  field  U(t,  x,  y,  z)  is 

V^U  - v"*  - 0 6(y)6(z)6(t),  U = 0,  t < 0.  (4) 

Here  is  the  three-dimensional  Laplacian  and  6 denotes  the  Dirac 
delta  function.  At  each  point,  x = ^,  y = z = 0,  on  the  source- 
receiver  line,  an  impulsive  source  is  set  off  and  the  backscattered 
field  at  (t,  5,  0,  0),  denoted  by  Uj(t,  5),  is  observed.  Each  such 
experiment  (i.e.,  for  each  c)  is  set  off  separately.  In  each  case, 
the  "time  clock  of  observation"  is  reset  so  that  t = 0 corresponds  to 
the  time  when  the  pulse  is  set  off. 

This  formulation  dispenses  with  the  upper  surface  as  a boundary. 
Alternatively,  one  could  replace  (4)  by  a homogeneous  equation  in 
z > 0 and  prescribe  a boundary  condition  at  z = 0 which  contains  the 
source  term.  In  the  simplest  model  of  this  type  one  would  prescribe 
9U/az  as  an  impulse  at  each  surface  point.  The  consequent  change 
induced  by  this  simplest  model  is  to  introduce  multipliers  of  2 and  4 
in  the  results  stated  below,  while  not  changing  the  substantive  results 
at  all.  Thus,  we  proceed  with  (4)  as  the  governing  equation. 


With  c denoting  a constant  reference  speed,  we  define  the 
variation  in  v"^(x,  z)  by  the  equation 

v''(x,  z)  = c‘"[l  + a(x.  z)]  . (5) 

By  using  the  methods  in  Bleistein  and  Cohen  (1977),  the  following 
integral  equation  may  be  derived  for  a(x,  z): 


(6) 


Here,  H is  the  Heaviside  function, 

p = /(x-^)^  + z^ 


(7) 


and 

T 

e(T,0  = - (4ttc)' J dt(T-t)U3(t,d.  (8) 

0 

The  integral  equation  can  be  solved  by  transform  techniques. 

The  result  is 


♦ 


See  Appendix  A. 
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•U^(t,C)exp{2ik^(x-C)  - 2ik^z  + iairl  ; (9) 


This  is  a formula  which  reproduces  a(x,  z)  by  direct  processing  of 
the  observed  data  itself.  In  fact,  this  is  not  the  formula  to  be 
used  for  computer  implementation.  Instead,  the  integrand  is 
multiplied  by  -2ik  to  yield  a formula  for  a'  = 9a/8z. 

3 

For  a layered  medium,  a(x,  z)  is  a constant  in  each  layer 
while  cx'  is  proportional  to  a sum  of  Dirac  delta  functions  which 
peak  on  the  interfaces  between  the  layers.  Since  band  limited 
delta  functions  are  much  easier  to  recognize  than  band  limited 
step  functions,  it  is  far  more  desirable  to  process  data  for 
n'  than  for  a.  An  analytical  development  of  these  ideas  can  be 
found  in  Mager  and  Bleistein  (1976).  In  particular,  it  is  shown 
in  that  paper  how  to  estimate  the  magnitude  of  the  jump  in  a at 
a discontinuity  in  terms  of  the  output  of  the  band  limited  deri- 
vative, a' . It  is  this  analysis  that  provides  our  velocity  esti- 
mates from  this  band  limited  data. 


a{x,  z)  = ^ jdE,  j dk_  j dk^  j 


Examples 

To  take  account  of  the  constraints  cited  earlier  in  this  paper, 
the  multiple  integration  in  (9)  is  carried  out  as  follows.  The 
time  domain  is  truncated  at  the  maximum  time  of  "reliable"  observa- 
tions UgCt,  c).  The  step  size  in  the  time  interval  is  the  sampling 
rate.  The  integrals  in  wave  vector  are  restricted  to  an  annulus 
consistent  with  the  frequency  constraint,  i.e., 

2TTf  Zirf^ 

^ < /in=  + k 2 < -r—  . 

C 1 3 C 

The  integration  over  ^ is  truncated  to  the  interval  (-L,  L)  and 
the  step  size  AC  (in  this  paper,  100  ft.)  is  used  for  this 
quadrature. 

A computer  algorithm  has  been  developed  by  the  authors  to 
generate  a'  in  accordance  with  (9)  and  the  discussion  above. 

Synthetic  data  U5(t,  C)  was  generated  for  various  subsurface  profiles 
and  then  the  data  was  processed  to  produce  a'. 

In  Figures  2,  3,  and  4 are  the  results  for  profiles  of  the  type 
in  Figure  1 with  0 being  15°,  30°,  and  45°  respectively.  The  angle  of 
inclination  of  the  interface  is  accurately  reproduced  in  the  region 
bounded  by  the  dashed  lines  of  Figure  1.  The  velocity  c in  this  case 
was  5,000  ft. /sec.,  Ac  was  250  ft. /sec.  A typical  value  of  Ac 
estimated  from  the  output  was  249.7  ft. /sec.  In  this  type  of 
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configuration,  the  error  is  virtually  linear  in  Ac  and  hence  these 
percentage  errors  will  be  maintained  for  much  highe"  values  of  Ac. 

For  these  examples,  one  can  generate  the  synthetic  data 
analytically  and  calculate  the  integral  in  (9)  asymptotically,  under 
the  constraint  (11).  The  result  of  that  analysis  is  completely  con- 
sistent with  the  numerical  output,  both  as  regards  the  location  of 
the  interface  and  the  estimate  of  Ac. 

Figure  5 shows  a time  section  for  an  anticline.  Figure  6 is  the 
result  of  direct  inversion  and  gives  virtually  the  exact  model, 
which  was  a circular  anticline.  For  example,  the  top  of  the  actual 
circle  was  located  at  1,500  ft.  and  our  output  yielded  an  estimate 
of  1,502  ft.  The  actual  jump  Ac  at  this  point  was  again  250  ft. /sec.; 
the  estimate  from  our  output  was  247.1  ft. /sec.  In  the  middle  of 
the  flanking  plane  (x  = 2,700  ft.),  the  actual  location  and  jump 
were  2,000  ft.  and  250  ft. /sec.,  respectively;  our  estimates  were 
2,000  ft.  and  248.1  ft. /sec. 

Figure  7 is  the  time  section  for  a circular  syncline.  Figure  8 
is  the  result  of  our  direct  inversion  procedure.  Depth  and  Ac  at 
the  bottom  of  the  circle  were  3,472  ft.  and  250  ft. /sec.;  direct 
inversion  yielded  3,472  ft.  and  249.3  ft. /sec.,  respectively.  At 
X = 3,600  ft.  (on  the  flank  plane)  the  depth  and  ac  values  were 
2,000  ft.  and  250  ft. /sec.;  direct  inversion  yielded  2,010  ft.  and 
253.6  ft. /sec. 
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Figure  9 shows  a subsurface  configuration  for  which  synthetic 
data  was  graciously  provided  to  us  by  the  research  group  at  Marathon 
Oil.  The  time  section  provided  by  them  is  depicted  in  Figure  10; 
Figure  11  is  the  result  of  our  direct  inversion  procedure  and 
Figure  12  shows  our  estimate  from  the  output  of  various  relevant 
quantities  in  the  model.  The  lower  two  sets  of  numbers  exhibit 
errors  of  less  than  4%,  while  above  that  level,  the  errors  are  less 
than  1%. 

Automatic  gain  control  was  applied  to  the  input  and  the  data 
was  again  processed  by  our  inversion  procedure.  The  output  was 
indistinguishable  from  Figure  11.  In  this  case,  the  output  does 
not  provide  velocity  estimates,  and  thus  our  inversion  procedure 
becomes  a migration  procedure. 
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Conclusions 


We  have  derived  an  approximate  solution  to  the  inverse  problem 
for  the  velocity  in  an  inhomogeneous  medium  which  supports  acoustic 
waves.  The  approximations  made  are  often  used  in  modeling  the 
inverse  problem  in  seismic  exploration  and  can  be  found  also  in  the 
references  cited  earlier,  as  well  as  in  many  other  papers.  A com- 
puter implementation  on  synthetic  data  under  a number  of  realistic 
constraints  has  also  been  carried  out. 
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APPENDIX 


The  derivations  of  equation  (6)  will  be  presented  here.  To 
begin,  the  function 

Uo(t,  X,  y,  z)  = 6(t  - r/c)/(4TTr), 

r = ✓ x2  + yz  + zz  (Al) 

is  introduced.  This  is  the  Green's  function  for  the  "wave  operator" 
in  (4)  when  v is  replaced  by  c,  defined  in  (5). 

We  form  the  expression, 

U.(t  - t.  X -c.  y,  2)  [7*  - C-"  ||r]  U(t.  X,  y,  2) 

- U(t,  X.  y.  2)  [7“-  c-“  gr]  U.(t  - t,  X • C.  y,  2), 

and  integrate  over  space-time.  The  integral  can  be  shown  to  be 
equal  to  zero  by  Green's  theorem.  An  alternative  result  is  obtained 
by  noting  the  following: 

['7^  - I^J  Uo(t  - t,  X - j;,  y,  z)  = -6(x  - c)6(y)6(2)  6(T-t). 

Uo  = Uo^  = 0,  t = t;  (A2) 
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fF”]  ' c)'5(y)6(z)6(t)  + Jr  Ip  . 

U = U^  = 0.  t = 0.  (A3) 


Thus,  it  follows  that 


00  00 


•00  —00 


t,  X - 5,  y,  z)  |^-6(x  - 5)6(y)6(z)6(t) 


+ 


a 3^U  I 

'W  J 


+ U(t,  X,  y,  z)6(x  •c)6(y)6(z)6(T  - t)|  . (A4) 

By  using  the  Dirac  delta  functions  and  the  initial  and  final  condi- 
tions in  (A2,  3)-thts  result  may  be  simplified  to  the  following: 


dxdydz  a 


(A5) 


= - c2U(t,  C,  0,  0)  = - U3(t,  0-  (A5) 

Here,  the  last  expression  defines  Uj(t,C)  as  used  in  this  paper. 
Integration  by  parts  in  (A5)  and  use  again  of  the  initial  and  final 
conditions  in  (A2,  3)  places  the  second  time  derivative  aV^t^ 
on  Uo  instead  of  on  U.  This  can  then  be  replaced  by  a second  time 
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derivative  with  respect  to  t : A double  integration  with 

respect  to  t and  an  interchange  of  orders  of  integration  then 
yields 


oo 


/./ 


dxdydz  a Uo(t  - t,  x - C.  y.  z)U(t,  x,  y,  z) 


c*  (t-  t)U5{t,  C)  dt. 


(A6) 


If  (A3)  is  formally  solved  by  perturbation  methods,  then  to  leading 
order 


U(t,  X,  y,  z)  = Uo(t,  X - y,  z) 


(A7) 


with  the  correction  term  of  order  cx.  Substituting  (A7)  yields  the 
equation, 

T oo 

y*dt ^ dxdydz  a(x,  z)  Uo(t-  t,  x -C.  y,  z)Uo(t,  x - C,  y.  z) 


0 -00 


I 

= -c^y*  (t-  t)U5(t,  £)dc  , 


(A8) 


with  the  correction  to  the  integrand  on  the  left  of  the  order  of  a^, 
which  is  negligible  for  a small.  With  U given  by  (Al),  the  t and  y 
integrations  can  now  be  carried  out  to  yield  (6). 
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FIGURE  10 


FIGURE  11 


